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ABSTRACT 

A two-dimensional finite-difference code to solve the 
BGK-Boltzmann equation has been developed. The so- 
lution procedure consists of three steps: (1). Transform- 
ing the BGK-Boltzmann equation into two simultaneous 
partial differential equations by taking moments of the 
distribution function with respect to the molecular ve- 
locity u 2 , with weighting factors 1 and it?. (2). Solving 
the transformed equations in the physical space based on 
the time-marching technique and the four-stage Runge- 
Kutta time integration, for a given discrete-ordinate. 
The Roe's second-order upwind difference scheme is used 
to discretize the convective terms and the collision terms 
are treated as source terms. (3). Using the newly calcu- 
lated distribution functions at each point in the physical 
space to calculate the macroscopic flow parameters by 
the Modified Gaussian quadrature formula. Repeating 
Steps 2 and 3, the time-marching procedure stops when 
the convergent criteria is reached. A low-density nozzle 
flow field has been calculated by this newly developed 
code. The BGK Boltzmann solution and the experimen- 
tal data show excellent agreement. It demonstrated that 
numerical solutions of the BGK-Boltzmann equation are 
ready to be experimentally validated. 

INTRODUCTION 

Low density nozzles are often used as a propulsive 
system for the trajectory control of satellites and other 
spacecrafts in orbit. One of the distinct characteristics of 
the low-density nozzle flow is that the expanding gas is 
highly rarefied with a thick boundary layer near the wall 
such that the distinction between the boundary layer and 


t Assistant Professor. Member AIAA 
+ Professor, Member AIAA 
* AST, Senior Member AIAA 

Copyright ©American Institute of Aeronau- 
tics and Astronautics. Inc., 1995. All rights 
reserved. 


the core flow disappears. Thus the fluid dynamics in 
a low-density nozzle falls into the transitional regime. It 
has been well known that the conventional Navier-Stokes 
equations fail to predict flow characteristics accurately in 
the transitional flow regime^ 1,2 ]. Usage of improper con- 
tinuum methods will lead to significant erroneous esti- 
mates of the nozzle flow characteristics^. A rigorous nu- 
merical treatment of the flow field within the low-density 
nozzle requires the solution of the Boltzmann equation. 

The Boltzmann equation is a highly non-linear 
integro-differential equation to describe flow characteris- 
tics in the molecular level with the molecular distribution 
function as the only dependent variable. Direct numeri- 
cal solutions of the Boltzmann equation for fluid flows in- 
volves the calculation of the molecular distribution func- 
tion as well as the collision integral at each velocity point 
in a three-dimensional infinite velocity spaced. 

Many research efforts have been reported in di- 
rectly solving the non-linear Boltzmann equation using 
finite-difference methods. The Monte-Carlo-Integration 
technique (MCI)^ 5,6,7 ] has been exclusively used to eval- 
uate the collision integrals. The shortcoming of the 
MCI method is the large consumption of computer time, 
and the MCI procedure must be executed at each time 
step. Without loss of generality, the linearized BGK- 
Boltzmann equation has been solved successfully with 
large Knudsen number^ 8,9 ]. The advantage of solving 
the BGK-Boltzmann equation is that it requires much 
less computational efforts and computer time, and still 
preserves the flow characteristics of the original Boltz- 
mann equation, as long as the Knudsen number is rela- 
tively large. For engineering applications, the solutions of 
the macroscopic parameters are required. These parame- 
ters are derived from the molecular distribution function. 
The discrete-ordinate method^ 8 * 9 ] is based on the con- 
cept that the integration over an infinite velocity space 
can be reduced to integration over a finite number of dis- 
crete points by selecting a appropriate integration for- 
mula. The transformed finite difference equations from 
the Boltzmann equation can be quantitized in the veloc- 
ity space, and the macroscopic properties can be obtained 
though this quantitized velocity space. 

The objective of this research is to develop a new 
finite-difference method to solve the BGK-Boltzmann 
equation for nozzle flows in the transitional flow regime. 
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The new methodology combines the Discrete-Ordinate 
(FDDO) technique^ 9 ! with a second-order upwind finite- 
difference method^ 10 ! and the four-stage Runge-Kutta 
time integration scheme. This paper describes the nu- 
merical procedure and presents solutions of low-density 
nozzle flows. 

MATHEMATICAL FORMULATION 

A two-dimensional BGK-Boltzmann equation, with- 
out external force, is solved in a Cartesian coordinate 
system 

% + Ux % + u »% = Ac{f °- f) (1) 

where f 0 is the Maxwellian distribution function, A c is 
the collision frequency, and 
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m is the molecular weight, n is the molecular number 
density, u x ,u y are the molecular velocity components in 
the x— and y— directions, respectively. It is assumed 
that the viscosity coefficient p and molecular mean free 
path Xi are 
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The subscript 0 indicates the reference state, and the 
exponent w depends on the flow media. Equation (1) is 
simplified as 


9/ 9/ 9/ _ 16 1 n ( RT 

dt Ux dx Uy dy 5 Xi n 0 \ 2 tt 


(2 jrRTy 


- Exp 


t- 


V 2 RT 




( 6 ) 


In order to reduce the number of independent variables, 
the BGK Boltzmann equation is integrated with respect 
to u, , with weighting factors 1 and u 2 z . Equation (1) 
becomes 
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Taking a one-one correspondence transformation in the 
velocity space as well as in the physical space 


u x = V sin d) 
u y = V cos 0 

0 = tan 


-1 /«£ 
\ Uy 


equations (7) and (8) become 


where 
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with Jacobian J (x, y;£, rj) and coordinate transformation 
coefficients x^,x n ,y^,y n . 

NUMERICAL PROCEDURE 

The macroscopic flow parameters can be obtained 
through the integration of the distribution function /. 
For example, the number density, A ’[x.y.t) in two- 
dimensional cartesian coordinates can be obtained from 
the zero-moment equation, 

N{x,y,t) = 

f°° f°° f°° ri (20) 
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where (u r ,i/y,u c ) is the molecular velocity vector in the 
velocity space. 

Substituting equation (9) into equation (20). it be- 
comes 

roo r OC 


/*OG /*OC 

N(x,y,t) = / / g(x,yA.u I .u,j)du r du y (21 
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where 


Using polar coordinates for velocity space to rewrite 
equation (21). i.e.. 



u T — 1' sin 6 


( 22 ) 

u y = V cos 0 (23) 

0 < I'm. 0 < 0 <2z 


Then, equation (21) is transformed into 


N(x,y,t) 



V sin 4>g{x, y,t. V', 4>)dodV 


(24) 


Since the molecular distribution function, /, is in an ex- 
ponential form in nature, it is reasonable to assume that 
g(x,y,t) is a function of the exponential type. Then, 
equations (24) can be rewritten as 


N(x,y, 
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\\\ = 0.379530780467479£ - 02 
W 2 = 0.213680829692996F - 01 
W 3 = 0.559585707892949 £ - 01 
W 4 = 0. 958716826650700 £ - 01 
W 5 = 0.116908207001337£ + 00 
W 6 = 0. 102936301 287559£ + 00 
W 7 = 0. 646824672793000 £ - 01 
W 6 = 0.283191 162204620£ - 01 
W 9 = 0.836264802590032£ - 02 
W 7 io = 0.159773621 113803£ - 02 
W n = 0.1870134659 16242£ - 03 
W 12 = 0. 12439355056 1664£ - 04 
W 13 = 0 .420846696 1 87 155 £ - 06 
W u = 0. 605184708943963 £ - 08 
W 1S = 0.264340659193899£ - 10 
W 16 = 0.15245941 1718563£ - 13 


where 
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and 

W k = W k e v ?,k = 1,2, 3, - . 16 
The root \\ is called the Discrete-Ordinate. 


Equation (25) can be integrated effectively using 
the modified Gauss-Hermite quadrature formula. This 
quadrature formula replaces the semi-infinite integral 
by a finite number of integration points, called “roots 
Vfc, k = 1,2, • • n”, and transforms equation (25) into 

n 

N(x, y, t) = '£w k P(V k ) (27) 

k-l 


where W k is the corresponding weighting factor for root 
V\. For lj = 1, n = 16, the root and its weighting factors 
are given as 


Vi = 0.477579953723861 £ - 01 
V 2 = 0.1 57564360925804 £ + 00 
V 3 = 0.323655656470272£ + 00 
V 4 = 0.539 147354 111002£ + 00 
V 5 = 0.797005397275377£ + 00 
Ve = 0.109095830650419£ + 01 
W = 0.141597596974798£ + 01 
V 6 = 0.176843702942131£ + 01 
V 9 = 0. 2 146 1499609 1 144 £ + 01 
Ho = 0.254836565149444£ + 01 
V'n = 0.297589659136340£ + 01 
V'i 2 = 0 .343 1483867 15786 £ + 01 
f 13 = 0.392069411852247£ + 01 
f 14 = 0.4454 12057238520£ + 01 
V is = 0.505367426854 191 JE7 + 01 
V 16 = 0. 577847884687290 £ + 01 


For any selected discrete-ordinate V k , 
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In this study, a 16 point discrete-ordinate method is 
adopted, N = 16. 

Introduce the following dimensionless variables 
x y n 

x = n = — 

a a n 0 

u x u y - T 

Ur = T 7 " ’ = T = — 

Vo To 

4 = M t= 1 

Vo d/-/2RTo 
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G = 


GV o 2 
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where d is the characteristic length of the flow field, the 
subscript 'O' represents the reference state, and 

Vo = yj 2 RTq , 


then, equations (28,29) become 



(34) 

(35) 


A successful numerical scheme for solving equations (34) 
and (35) should be fast and reasonably accurate. An 
assessment of numerical techniques for solving these 
equations concludes that the explicit Runge-Kutta four- 
stage time-integration scheme with second order upwind 
flux differencing^ 10 ! performed very well with regard to 
handling the nonlinearity and the convergence speed. 
Rewrite equations (34) and (35) into 
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and the Minmod operator is defined by 


( 38 ) 


where Q represents g k and h k in above equations, E is 
BkQ, F is CjbQ, 5 represents source terms. The four- 
stage Runge-Kutta integration technique is applied to 
solve equations (34) and (35). The solution procedure 
can be summarized as follows. 

For a given discrete-ordinate V ky equations (34) and 
(35) are solved, with 0 < <j> < 27 r, for every points ' 
in the physical space using the designated Runge-Kutta 
scheme. A characteristic-based spatial discretization us- 
ing second-order upwind flux-differencing is adopted to 
discretize the convective terms. The explicit spatial op- 
erator may be written in terms of numerical flux vectors 
treated on a locally one-dimensional basis. The residue 
R(Qij) can be written as 



The superscript (2) represents the second-order terms. 
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where 


Minmod[x,y] = 5pn(x)Max (0, |xj. y s£n(x)}) 

(39) 

The values of c is chosen based on the type of schemes, 
for second-order upwind scheme, c = —1.0. The constant 
ft is a compression parameter which is restricted to the 
range 1 < 3 < (3 - e)/(l - c) with (3 = 6 when c = 1. 

This process repeats for the entire discrete- 
ordinates. The nondimensional macroscopic parameters, 
such as velocities, temperature, and number density, are 
updated using the new values of g and h. at each point 
in the physical space, based on the Modified Gaussian 
quadrature 

16 y 2ir 

W k g k do 

k= 1 *'° 

2ir 

W k V k sin og> : do 

- 1 A [ 2 * 

U y = - > / IVjtU cos og k do 

n k=i Jo 

F = 3 7>Ej 0 W k(h k + l* 2 g k )do - n (r* + F 2 ) 

The integration over 0 is obtained by the Simpson's rule. 
Currently, 16 points are used in o direction. 

Three types of boundary conditions are implemented 
for the benchmark solution. 

(1). For a given discrete-ordinate l>. Reservoir Bound- 
ary Condition: 
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9k = 



hk = 



where ~ represents the non-dimensional quantity. 

(2). Symmetric Boundary Condition: 

Along the nozzle and reservoir centerline, symmetric con- 
dition is applied. For a given discrete-ordinate V*, 


9k\y = o(* “ <t>) = 9k\y = o{<t>)> 

0<^<f 


7T 

hfc|y=o(^ “ hk ly > 

0<«<j 

Qk Jy = 0 (3?T 0) -- 9k |y = o(^)> 

. 3tt 

7T < <f> < 

“ “ 2 


3 7T 

^jfc|y=o(37r — 0) “ hk |y = o(^)) 
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(3). Wall Boundary Condition: 

In order to specify the interaction between molecules and 
the wall surface, diffuse reflection is assumed, i.e. for 
a given wall temperature T w , and discrete-ordinate V*, 
molecules which strike the surface are then emitted ac- 
cordinf to a Maxwellian distribution, 
n - 

9k = —£~ e T “ . for (14 • n) > 0 

7T lu) 


hk — 2 9 k ) 


for (V* • n) > 0 


where, n is the inward normal vector at the wall surface. 
Since the wall number flux h w is unknown a priori, it is 
determined by the previous time-step value, 

h w = -2(4^)* [ [ (V-n)gVd<t>dV, for (V n) < 0. 

T w Jo J<j> 


Repeating above procedures, the time-marching proce- 
dure stops when convergent criteria are satisfied for all 
the macroscopic parameters. Currently, if the summa- 
tion of relative error of density, velocity and temperature 
over the entire physical space is less than 10“ 5 , the solu- 
tions are considered “converged”. 

RESULTS AND DISCUSSION 

Numerical results of low-density nozzle flow was 
obtained with the developed BGK-Bo!tzmann-2D code. 
The nozzle is placed between two reservoirs which repre- 
sent the upstream boundary condition and downstream 
boundary condition near the nozzle exit. The upstream 
reservoir pressure and temperature are 0.01 psi and 900 
R. The Knudsen number of the Argon gas in the cham- 
ber is 0.1, while the Knudsen number increases to 20 at 
the nozzle exit. The downstream reservoir pressure is 
0.0001 psi. Figure 1 shows the grid setup. Block logic 
is implemented in the code so that it has the ability to 
calculate flows with complex geometry. A more sophisti- 
cated grid network is under development and results will 
be addressed in the conference. The pressure ratio be- 
tween the two reservoirs is selected to be 100:1 for this 


case. The upstream and downstream reservoir tempera- 
ture is set at 900 R. The isothermal nozzle wall is set at 
temperature of 540 R. The Knudsen number, defined by 
the ratio of upstream reservoir mean free path to nozzle 
inlet diameter is selected to be 0.1. Figure 2 shows the 
velocity vectors inside the nozzle. Figure 3 shows the 
Mach number contours inside the nozzle. These conver- 
gent solutions were obtained with 8 hours of CPU time 
on the CRAY YMP. 

More extensive numerical simulations will be per- 
formed for the final presentation in the conference. As a 
result, a refined analysis of low-density nozzle flow char- 
acteristics will be summarized in the conference. 
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Velocity Vectors inside the Nozzle 
Reservor Pressure Ratio: 100:1 
Knudesn Number: 0.1 
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Figure 2. Velocity vectors inside the nozzle. 
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Mach Number Contour inside the Nozzle 
Reservor Pressure Ratio 100:1 
Knudesn Number: 0.1 


Figure 1. Crid setup for flow between two reservoirs. 

Scale is non-dimensionalized with respect to 
nozzle throat diameter. 
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figure 3 Mach number contours inside the nozzl 







